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1.   Introduction 

With  the  ongoing  development  of  larger  and  faster  computers,  global  predic- 
tion using  primitive  equation  models  incorporating  complex  physics  has,  since 
1975,  become  feasible  (i.e.  the  European  Mid  Ranii.e  Model,  NMC's  spectral 
model,  the  U.S.  Navy  NOGAPS  model). 

Based  on  work,  by  Rossby  dating  as  far  back  as  1938  and  the  subsequent  work 
of  Cahn  (1945)  concerning  geostrophic  adjustment,  it  r*as  realized  already  in 
the  early  1950' s  (Hinkleman,  1951),  that  initializing  a  primitive  equation 
model  with  some  form  of  a  quasi-geostrophic  balancec  v:ind  (in  Hinkleman' s 
early  work  -  just  the  geostrophic  wind  itself)  would  be  adequate  to  obtain 
useful  forecasts  since  the  geostrophic  adjustment  mechanism  permitted  by  the 
primitive  system  would  be  able  to  keep  the  inert lal-gtavity  wave  amplitudes 
small  relative  to  the  Rossby  wave  amplitudes. 

Due  to  increased  understanding  of  how  to  properly  control  the  development 
of  nonlinear  instability,  practical  primitive  equation  models  were  developed 
in  the  1960's  first  for  general  circulation  study  (to:  which  there  was  really 
no  initialization  problem)  and  then  for  operational  .'Jorecasting .   Increasingly 
refined  methods  for  obtaining  a  balanced  initial  Stat*?  have  been  developed. 
Recent  among  them  have  been  the  variational  methods  iad  the  normal  mode 
method.   In  the  latter,  the  inertial- gravity  waves  and  the  mass  and  wind 
tendencies  associated  with  them  are  explicitly  removed  from  the  initial  fields 
(Machenhauer ,  1977). 

In  the  late  1960 's  another  and  quite  simple  initialization  method,  now 
known  as  dynamic  initialization  was  proposed.   Reference  for  it  can  be  found 
in  Nitta  and  Hovermale  (1969),  Miyakoda  and  Moyer  (1968)  or  Winninghoff  (1968, 
1973a,  1973b).   It  enjoyed  some  attention,  but  because  of  its  need  for  a  large 
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amount  of  computer  time,  interest  in  it  faded  especially  in  comparison  with 
the  more  efficient  normal  mode  technique. 

The  purpose  of  the  work  which  is  reported  on  in  the  present  paper  was 
originally  to  extend  the  work  in  1973b.   Motivation  also  comes  from  the  more 
advanced  state  of  computer  technology  at  present. 

In  1973b,  the  author  used  a  free  surface  global  model  with  a  10  mesh  size 
to  investigate  a  dynamical  initialization  procedure.  Results  clearly  indicated 
that  restoration  to  the  mass  field  in  high  latitudes,  to  the  wind  field  in  low 
latitudes  and  a  blend  of  the  two  between  30   and  60  was  very  satisfactory  for 
the  purpose  of  reducing  the  amplitudes  of  the  so  called  "noise"  waves.   Then 
an  investigation  was  made  concerning  incorporation  of  heating  and  friction 
into  the  procedure.   Here  rather  idealized  circular  symmetric  tropical  pertur- 
bations were  used.   The  result  clearly  indicated  the  feasibility  of  doing  this 
when  the  equations  were  run  forward.   The  author  was  using  a  forward- backward 
iterative  procedure.   During  the  backward  part,  however,  whether  restoration 
be  done  to  the  total  wind  or  just  the  nondivergent  part  was  not  clearly  deter- 
mined.  The  model's  mesh  size  was  too  crude  since  the  computer  used  at  that 
time  was  very  small  compared  to  those  available  today. 

The  use  of  global  prediction  models  must  lead  to  increased  interest  in 
some  of  those  previously  unwanted  inertial-gravity  type  waves  especially  as 
they  affect  mesoscale  weather  patterns  in  the  midlatitudes  and  as  they  are 
present  in  the  form  of  equatorial  waves  of  mixed  Rossby-gravity  type.  Forcing 
by  heating  may  of  course  be  involved  for  all  these  types  of  waves  so  there- 
fore even  the  normal  mode  method  of  initialization  may  be  limited  and  even 
more  so  near  the  equator,  where  the  separation  of  the  waves  based  on  frequency 
is  not  feasible. 
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In  Che  second  part  of  the  earlier  work  (1973b),  the  author  did  begin  to 
consider  the  inclusion  of  mass  forcing  and  friction  into  the  dynamic  initiali- 
zation procedure  using  a  global  free  surface  model  but  the  model  was  not  only 
Loo  crude,  having  only  the  10  degree  mesh  length  but  also  it  did  not  have  the 
best  spatial  staggering  of  points  for  the  purpose  of  simulating  geostrophic 
adjustment . 


II.   Model  Experiments 

The  model  used  in  this  work,  is  still  the  free  surface  primitive  equation 
model.   It  describes  an  incompressible,  homogeneous  fluid  with  a  flat  lower 
boundary,  a  free  upper  surface  and  is  subject  to  rotation.   It  is  applied  on  a 
sphere.   The  finite  difference  equations  are  based  on  the  energy  conserving 
requirement  proposed  by  Lorenz  and  applied  and  developed  by  Arakawa.   They  are  ; 
written  in  momentum  flux  form  thus, 

i«<S>+{r  <£■>+*?  <£'>-&  + 

at  Km'      ax  VN      '      scp   M     '      1mij 

/31  31N1l  d)9d)„  /on 

3A   N  3cp  M  N    3(£  f 

It  4> +  ix  <r> + k <f>  ■&  (3) 

where  1/M  =  a  cos  cc;  l/N  =  a  and  a  =  radius  of  the  earth,  co  =  latitude,  X  e 
longitude,  $  =  geopotential ,  u  =  zonal  component  of  the  wind,  v  =  meridional 
component  of  the  wind,  f  =   coriolis  parameter,  Kf  e  is  a  Newtonian  friction 
coefficient  andjcr'  represents  seme  as  of  yet  unspecified  source  or  sink  of 
mass  (a  simulated  "heating"  or  "cooling"). 

The  finite  difference  equations  are  now  written  using  the  Arakawa  scheme  C 
(Winninghoff ,  1968;  Arakawa  and  Lamb,  1977)  which  is  superior  for  describing 
geostrophic  adjustment.   For  these  experiments  a  72  degree  wide  sector  of  the 
sphere  is  used  and  the  increments  of  longitude  and  latitude  are  4  degrees.   A 
six  minute  time  step  is  used.   Therefore,  because  of  the  convergence  of  the 
meridions  a  Fourier  smoothing  of  the  zonal  pressure  gradient,  the  zonal  con- 
vergence in  equation  3,  and  the  zonal  components  of  the  horizontal  advection 
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in  Eqs  (1)  and  (2)  is  used  (Arakawa,  1971),  to  satisfy  the  linear  stability 
criterion.   Though  the  leapfrog  time  step  is  normally  used  during  the  regular 
forecast  runs,  an  Euler  backward  step  is  inserted  every  5  time  steps  mainly  to 
control  any  solution  separation  caused  by  use  of  mass  source  terms.   During 
the  initialization  phase,  to  be  explained  later,  only  the  Euler  backward  steps 
are  used. 

To  review,  during  the  dynamic  initialization  phase  the  model  is  integrated 
alternatively  forward  and  backward  for  a  few  time  steps  until  the  fields  of 
mass  and  wind  appear  to  settle  into  near  balance.   It  is  during  :hi  s  procedure 
that  the  Euler  backward  scheme  is  used  exclusively  since  it  has  damping 
effects  on  the  high  frequency  motions.   In  addition,  a  constraint  or   the  ::orm 


n+1  .    n+1   „  ,      n+1    o, 
'"a 


a  =  a+     -  K_  (a*     -  a  )  (4) 


can  be  applied  to  the  u,  v,  or  $    fields  or  all  of  them  where  the  K  's  might  be 

latitude  dependent.   (A  variation  could  be  to  have  the  K  's  scale  impendent 

a 

explicitly.   In  the  east-west  direction  Fourier  analyses  could  be  :sed  while 

in  the  north-south  direction  some  other  technique  would  be  needed-)   Where  K 

A   n+1     n+1    ,  ,     _    .    n+1     o  .  .  ,  .   . ,   .  .  .  ,    ,      _ 
=  0,  a  =  a*    an<i  where  K  =  1,  a    =2  whicn  is  the  inxtial  value.   To 

give  full  flexibilitv  K  's  are  usually  limited  to  a  value  no  larger  than  about 

a 

.5.  This  whole  method  is  quite  simple  and  allows  the  equations,  in  whatever 
form  they  are  cast,  to  approach  a  consistent  semi- balanced  state  through  the 
internal  adjustment  process  itself. 

The  foundation  of  much  of  the  work  reported  on  here  is  the  theory  of  geo- 
strophic  adjustment.   From  that  theory,  an  important  parameter  is  obtained  - 
the  radius  of  deformation,  /$"/f,  which  is  the  important  length  scale  parameter 
separating  what  can  be  considered  large  or  small  scale  relative  to  the  adjust- 
ment process.   In  Figure  1  is  plotted  the  radius  of  deformation  using  i>  =  9800 


-7- 


m  /sec  in  the  solid  line.   If  a  typical  scale  of  a  wave  is  defined  as  L/2tt 
note  that  not  only  are  all  scales  of  motion  small  near  the  equator,  but  even 
at  60  degrees  waves  of  6000  kilometers  or  less  in  wave  length  are  not  large 
and  thus  the  winds  are  still  an  important  source  of  initial  data.   In  the 
model  used  here  with  a  mesh  size  of  4  degrees  AX  at  60  degrees  is  only  222 
kilometers.   So  even  a  18  AX  wave,  the  longest  permitted,  is  about  4000  kilo- 
meters in  length  and  is  still  small  scale.   In  the  earlier  work  a  10  degree 
mesh  was  used  around  the  entire  360  degree  circle  so  therefore,  larger  scale 
waves  were  possible  at  higher  latitudes.   The  results  there  reflect  this.   The 
focus  in  this  work  will  be  directed  more  to  the  small  scale  because  of  the 
interest  in  the  equatorial  area.   The  broken  line  in  Figure  1  shows  the  scale 
of  the  18  AX  wave  (longest  permitted)  in  these  calculations. 
Experiment  1 

Here  the  purpose  was  to  verify  and  generalize  the  results  given  in  the 
second  half  of  1973b  for  the  "diabatic"  type  of  initialization. 

To  begin,  a  control  run  was  set  up  and  run  84  hours.   From  the  data  gene- 
rated 48  hours  into  that  run,  other  runs  were  made  using  variations  of  dynamic 
initialization  so  that  comparisons  could  be  made  of  the  new  forecasts  to  the 
original  control  run.   Usually  the  geopotential  at  some  point  was  used  for 
this  purpose.   The  change  in  the  free  surface  height  is  somewhat  like  the  tem- 
perature fields  in  a  baroclinic  atmosphere  and  thus  related  to  the  vertical 
motions  as  well.   Here,  however,  only  the  one  external  mode  is  involved. 

The  initial  data  for  the  control  run  was  specified  using 

<t>°  =  g[(1000  +  480  cos<£)  +  (100  sinX;L)61  +  (150  sin  nX  x2)<52J        (5) 
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Figure  1.   Radius  of  information  as  a  function  of  latitude,  H  =  1000  meters 
C solid  line).   Approximate  scale  of  18  A  wave,  the  longest 
pemicted  [  dashed  line). 
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Figure  2.   Temporal  evolution  of  $  at  a  selected  point  at  2°S,  84  hour  fr 
run. 
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for  the  geopotential  where  $°  is  given  in  m  /sec  .   Here  °P  is  the  latitude,  A 
is  the  longitude,  y  is  global  wave  number  5(1  for  the  limited  72  degree  wide 
grid)  ,  Xi  anc*  Xo  are  also  latitude  angles  but  are  modified  so  that  x-i  runs 
from  0  to  ir  for  Cp  from  +  62  degrees  to  +  2  degrees  and  x2  runs  from  0  to  u  for 
Cp  from  +  62  degrees  to  +  30  degrees.   Outside  these  ranges  6.  and  62  are  0. 
The  initial  velocity  components  are  (in  m/sec) 


u  -  -  ^p  (r->  '  v  =  ax  cf  >  •  (6) 

o  o 

-4    -1 
where  f  =  10   sec   .   (5)  gives  a  geopotential  increasing  into  the  tropics, 

but  with  a  slight  dip  near  the  equator  and  allows  for  one  wave  in  the  middle 

latitudes.   This  run  is  initialized  deliberately  rather  poorly  in  order  to 

excite  inertial-gravity  waves  of  sigrificant  amplitude.   In  addition,  a  mass 

source  term  is  applied, 


J4     =  (0.5  -  |sincp|)4»2At(.4'10)  (7) 


where  j^  is  added  to  $  at  each  rime  step.  This  amounts  to  about  .7  meters 
near  the  equator  tapering  off  to  0  at  30  degrees  N/S,  and  to  a  loss  of  about  .7 
meters  near  the  poles.  In  addition,  z.  point  source  of  about  2  meters  each  time 
step  is  allowed  at  a  point  located  at:  2  degrees  South  near  the  middle  of  the 
grid. 


Friction  is  included  with  L  =  3  x  10   •   Figure  2  displays  the  value  of  <j> 
at  the  point  where  the  point  mass  source  is  applied  and  Figure  3  shows  <J>  at  a 
point  at  46  degrees  South.   Note  especially  the  phase  and  the  amplitude  of  the 
inertial-gravity  oscillations  after  hour  48.   Amplitude  seems  to  average  about)! 
10  meters  or  about  1  percent  of  the  total  depth  and  the  two  most  prominent 
frequencies  seem  to  be  about  one  oscillation  in  6  hours  and  one  oscillation  inj: 
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Figure  3.   Temporal  evoluation  of  $  at  a  selected  point  at  46  S,  84  hour 
free  run. 


20  hours  both  near  the  equator  and  in  the  mid-latitudes.  At  46  degrees  however 
the  amplitudes  were  reduced  to  near  5  meters. 

Some  physical  understanding  of  why  such  waves  are  generated  may  be  sought 
in  examination  of  the  initial  conditions  and  the  nature  of  the  mass  source 
which  is  a  point  mass  source  at  2  S.   A  harmonic  analysis  in  the  east-west 
direction  does  indicate  a  maxima  in  wave  numbers  1,  2  and  4  in  the  i>  field 
near  the  equator  well  into  the  forecast.   These  have  wavelengths  of  about 
8000,  4000  and  2000  kilometers  respectively.   If  C,  the  phase  speed  of  gravity 


waves  near  the  equator  is  /$  this  is  /15000* or  about  120  m/sec.   The  frequency 
P  can  be  computed  from  L/C,  so  for  wave  numbers  2  and  4,  P  is  approximately  9 
hours  and  4.5  hours,  respectively.   For  wave  number  1  it  would  be  about  18 
hours.   In  the  initial  conditions  the  second  term  in  Eq  (5)  will  give  a 
wavelength  of  about  64   in  the  north-south  direction  (or  -  7050  kms)  and  the 
point  source  term  will  excite  waves  of  all  lengths.   So,  to  give  an  example  at 
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hour  60,  at  2°S  where  the  mean  <t>  is  14996  m  /sec  ,  the  amplitude  of  wave 

2    2  2    2 

number  1  is  44  m  /sec  ,  of  wave  number  2  is  11.8  m  /sec  ,  of  wave  number  3  is 

2    2  2    2 

1.32  m  /sec  ,  of  wave  number  4,  2.38  m  /sec  and  even  of  wave  number  9  is  .06 

2    2 
m  /sec  which  is  insignificant. 

Wave  number  1  is  predominant  in  the  midlatitudes  as  well.   Of  course  this 

comes  from  the  third  term  in  Eq.  5.   Its  length  will  vary  froii  about  7000  kms 

at  30   to  less  than  4000  km  above  60  .   Here  the  computation  of  C  should  be 


+  An  +  f2/K,2 


from  +  /gH  +  f  /K,   since  f  is  significant.   Also  its  dimension  in  the  y  di- 
rection  should  be  considered.   Using  45  a  value  of  C  of  about  140  m/sec  is 
obtained.   So  here  P  ~  (8.5  x  10  /140)   '  or  16  hours.   Uncertainty  in  this 
computation  is  sufficient  to  make  it  also  close  enough  to  explain  :he  longer 
observed  period. 

Next,  consider  the  first  term  of  Eq  (5).   The  mean  b   at  58   do  is   show  some 

2    2 
indication  of  a  very  long  oscillation.   Its  value  rises  to  12895.4  m  /sec  at 

2    2 
24  hours,  and  then  falls  monotonically  to  12743.2  m  /sec"  at  ?2  hoars.   The  82  I 

hour  run  was  not  long  enough  to  define  the  period  completely  but  a  1/2  period 
of  about  50  hours  seems  indicated.   If  the  first  term  in  Eq  (5)  is  differen- 
tiated with  respect  to  cp  to  obtain  an  expression  for  a  mean  U  due  to  this  term 

* 
a  value  of  about  5  m/sec  is  obtained  (480  (sin  cp)/f  a)  .   :\gain  assuming  a  L 

O  A 


of  4000  kilometers,  P  =  4,000,000  (m)/5(m/sec)  or  800,000  sec  which  is  of  the 
order  of  100  hours.   This  is  the  synoptic  scale  time  variation. 

Unfortunately  the  predominance  of  the  shorter  wave  is  not  explained.   It 
results  probably  from  imbalances  generated  by  the  forcing,  wave  number  inter- 
actions, and  possibly  difficulties  generated  near  the  poles  (here  wavelengths 
would  be  very  short) .   The  initial  imbalance  near  the  equator  is  quite  large  - 


a  is  the  radius  of  earth  -6.3  x  10  m  and  f   is  10    sec 

o 
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Figure  5a 


46  s 


Figure  4.   36  hour  forecasc  run,  temporal  evolution  of  $  at  a  selected  point 

at  2  S,  initialized  with  s em i~geo atrophic  wind  (solid  line),  $  from 
last  36  hours  of  free  run  (dashed  line).  Figure  4a,  same  as  Figure 
4  for  point  at  46  S.   Initial  $    in  forecast  run  from  free  time. 

Figure  5.   Same  as  Figure  4  except  *  replaced  between  30  S  and  30  N  in  the 
free  run  at  hour  48  with  a  constant  value  of  15060  m  /sec"  to 
initialize  the  forecast  run.   Figure  5a,  same  as  Figure  5  for  point 
at  46°S. 


in  fact  the  initial  conditions  were  set  up  deliberately  to  generate  noticeable 
noise . 

Figures  4,  4a,  5,  and  5a  show  the  <J>'  s  evolution  :.n  time  for  a  36  hour  prog- 
nostic control  run  initialized  from  data  taken  ."rom  48  hours  into  the  original 
84  hour  run.   Figures  4  and  5  show  the  1> '  s  at  2  degrees  South,  Figures  4a  and 
5a  the  <j>'s  at  46  degrees  South.   Figures  4  and  4a  are  for  a  run  where  the 
winds  were  initialized  from  Eqs  (6)  and  where  the  <J> '  s  were  taken  from  the 
control  run  itself  at  48  hours.   Figures  5  and  5a  are  for  a  run  where  the  t  at 

t=0  values  for  the  forecast  run  between  30  degrees  North  and  South  were  set  at 

2  ,   2 
a  constant  value  of  15050  m  /sec".   The  winds  again  came  from  Eqs  (6).   For 
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comparison  in  the  broken  curves  the  <|>  values  from  the  last  36  hours  of  the 
control  run  itself  are  plotted.   Note  the  pattern  of  the  noise.   In  neither 
ran  does  the  phase  of  any  inertial  gravity  motion  bear  any  resemblance  to  that 
indicated  for  the  final  36  hours  in  Figures  2  and  3.   Figures  4  and  4a  are 
somewhat  smoother.   Note  in  Figures  4  and  5  the  large  amplitude  wave  at  hour 
27  1/2  and  in  Figures  4a  and  5a  the  strong  maxima  at  hour  17  1/2.   There  are 
no  such  oscillations  in  the  broker,  curve. 

Figures  6  and  6a  show  results  at  the  same  two  grid  points  where  dynamic 
initialization  was  carried  out  for  12  complete  forward  backward  cycles  and  the 
cycle  length  was  12  time  steps.   No  mass  source  or  friction  terms  were  used  in 
the  initialization  but  they  were  applied  during  the  subsequent  forecast.   The 

F  's  of  equation  4  were  set  at  .5  for  K   and  K   from  30  degrees  South  to  30 

a       n  u      v  ° 

c  egrees  North,  and  K,  was  set  at  0  over  this  same  range.   North  or  south  of 

these  latitudes  K  and  K  varied  linearly  in  latitude  from  .5  to  .1  at  the 
u      v 

poles.   This  reflects  the  fact  (see  Figure  1)  that  even  the  longest  waves  per- 
mitted are  still  somewhat  small  scale  in  the  mid-latitudes.   In  fact  much 
l:esting  of  various  values  and  distributions  of  the  K's  yielded  the  best 
results  for  the  values  above. 

Comparison  of  the  curves  in  Figures  6  and  6a  to  the  broken  curves  now 
shows  a  very  important  result.   To  review,  the  reader  is  reminded  that  the 
tropical  t's  were  set  to  constant  values  at  the  beginning  of  the  dynamic 
initialization.   Unlike  the  previous  results  here  the  "real"  inertial  gravity 
motions  are  to  some  extent  recovered.   This  is  most  striking  in  the  phase  of 
the  waves  though  the  amplitudes  are  not  bad  either.   Still  the  amplitudes  in 
the  initialized  run  are  definitely  higher  than  the  control  run. 

Next  a  considerable  degree  of  testing  of  the  incorporation  of  the  mass 
source  and  friction  terms  into  the  dynamical  initialization  was  carried  out 
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Figure  6.   36  hour  forecast  run,  temporal  evolution  of  $  at  a  selected  point 
at  2  S,  adiabatic  dynamic  initialization  (solid  line),  $  from  last 
36  hours  of  free  run  (dashed  line).   Figure  6a,  same  as  Figure  6 
for  point  at  46  S. 

Figure  7.   36  hour  forecast  run,  temporal  evolution  of  *  at  a  selected  point 

at  2  S,  "diabatic"  dynamic  initialization  (solid  line),  *  from  last 
36  hours  of  free  run  (dashed  line).   Figure  7a,  same  as  Figure  7 
for  point  at  46  S. 


much  along  the  lines  of  that  done  in  the  earlier  work  (1973b).   Here  it  was 
reestablished  that  the  simplest  procedure  gave  the  best  result.   The  mass 
source  and  friction  terms,  formulated  exactly  in  the  same  manner  as  in  the 
free  run  were  allowed  to  operate  only  during  the  forward  part  of  the  cycle. 
Wind  restoration  was  done  during  the  entire  cycle  and  the  total  wind  was  used 
throughout.   It  was  found  that  this  was  superior  to  attempting  to  segment  the 
wind  into  divergent  and  nondivergent  parts  and  restore  to  one  part  or  the 
other  etc.  (This  may  be  due  to  computational  difficulties  of  solving  for  such 
winds  near  the  equator.) 
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An  operational  problem  with  this  procedure  is,  of  course,  that  the  mass 
source  and  friction  terms  must  be  allowed  only  to  operate  on  the  tendencies  of 
mass  and  wind  but  not  to  allow  any  accumulation  or  depletion  of  mass  or  energy 
during  the  initialization  procedure.   Merely  requiring  that  zonal  means  stay 
constant  is  not  sufficient  (as  in  1973b)  as  that  cannot  control  the  effects  of 
the  point  source  term.   So  what  was  done  here  vas  tc  aggregate  the  tendencies 
at  each  point  due  to  the  mass  friction  terms  during  the  forward  part  of  the 
cycle  and  then  to  subtract  them  out  completely  from  u,  v,  and  <j>  at  the 
beginning  of  the  backward  part  of  the  cycle.   This  is  not  quite  the  same  as 
running  these  terms  backward  (formally  possible  in  this  somewhat  artificial 
situation)  as  this  would  make  no  sense  physically.   Still  this  area  may  need 
further  research  to  be  applicable  to  a  "real"  atmosphere.   At  any  rate,  the 
results  shown  in  Figures  7  and  7a  are  for  this  modified  or  diabatic  initiali- 
zation procedure.   They  show  some  improvement  over  :he  results  in  Figures  6 
and  6a  in  the  amplitudes.   The  phase  is  still  well  ::orecast  and  the  amplitude 
is  improved. 

Some  further  tests  were  done  discarding  during  ':he  initialization  some  or 
all  of  the  "real"  mass  values  at  points  in  the  mid-latitudes .   Recall  from 
Figure  1  where  it  is  suggested  that  even  up  to  the  poles  these  runs  remain  to 
some  extent  within  the  small  scale  range  of  motions  though   the  18  AA  wave  is 
close  to  the  radius  of  deformation.   The  results  of  these  additional  runs  do 
indicate  that  some  mass  observations  are  needed.   The  case  with  no  mass  obser-' 
vations  showed  considerable  deterioration  again  with  respect  to  the  amplitudes 
of  the  forecast  inertial-gravity  waves.   Inclusion  of  some  mass  observations 
then  moved  the  result  almost  linearly  back  toward  that  of  Figure  7. 

The  importance  of  these  results  lies  in  the  implication  they  have  for 
tropical  initialization  at  the  present  time  and  eventually  for  initialization 
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of  mesoscale  motions  at  all  latitudes.   Referring  to  the  work  of  Matsuno 
(1966)  as  reviewed  in  Haltiner  and  Williams,  pages  99-102,  the  dynamics  of 
some  potentially  important  equatorial  waves  are  reviewed,   That  these  may  have 
profound  effects  under  certain  situations  has  been  emphasized  in  Webster's 
studies.   Though  difficult  to  verify,  this  type  of  wave  may  have  been  respon- 
sible for  the  extremely  heavy  precipitation  throughout  central  and  southern 
California  in  January  1969.   Thus  in  order  to  hope  for  any  capability  beyond 
the  4  to  6  day  range  with  numerical  models  better  initially t ion  of  Kelvin 
type  and  mixed  Rossby-gravity  type  equatorial  waves  may  be  crucial 

Therefore  several  experiments  were  set  up  to  examine  the.  dynamical  ini- 
tialization procedure  with  respect  to  the  larger  scale  in  spree  (not  relative 
to  the  radius  of  deformation)  Kelvin  and  Rossby-gravity  typo  v;aves    First,  the 
model  was  run  84  hours.   The  initial  conditions  were  taken  as  a  sinple  zonal 
flow  with  a  maximum  wind  of  10  meters  per  second  at  the  equator,  no  meridional 
flow,  and  the  heights  and  winds  in  geostrophic  balance  via  the  linear  balance 
equation.   Here 

u  =  10  cos  cp(m/sec) 

v  =  0.0     (m/sec)  (8) 

2r10  ,0n   10*.     2   ,  2,        2, 
M      2a       a 

where  <j>   =  5400  x  9.81  m  /sec  ,  a  e  radius  of  the  earth  (6.373  x  10  m)  ,  2ti    = 
14.548  x  10   sec   ,  cp  is  latitude  in  radians. 

A  Kelvin  wave  is  generated  by  means  of  a  simulated  mass  source  term.   It 
is  of  the  form, 

_7T_ 

H  =  A  sin  (nX)  cos(Y  Y)  (9) 


7 
where   A  =   25   m   /sec,    u   =   2-ir/L    (where   L   is    set    for  wave   number   1    for   the 
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Figure  8.   Free  run  with  forced  Kelvin  wave,  temporal  evolution  of  *,  U,  V  at 
selected  points  near  each  other  in  C  scheme  at  2  S  (V  on  equator) . 


model's  grid  which  would  correspond  to  wave  number  5  for  the  globe),  and  Y  is  a 
linear  function  ranging  from  -1  to  1  between  30  degrees  South  and  30  degrees 
North.   H  is  set  equal  to  0  in  the  mid-latitudes.   Friction  is  allowed  using 
the  same  formulation  as  earlier.   In  this  case  the  mass  source  is  stationary 
and  is  imposed  on  a  zonal  current. 

In  the  remaining  figures  (8  through  25a, b,c)  what  is  shown  is  first  a  plot 
with  time  at  the  point  at  2  degrees  South  of  the  u,  v,  and  $  values  and  then 
maps  of  the  u,  v,  and  <j>  fields  with  some  analyses  made  near  the  equator. 

For  the  first  84  hour  run  which  generates  a  Kelvin  wave,  Figure  8  displays 
the  temporal  plot  of  u,  v  and  <J)  and  Figures  9a,b,c  show  the  field  maps  at  12 
hours  into  the  run.   Figures  10a, b,c  then  show  the  field  maps  36  hours  later 
at  48  hours. 
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Next  an  initialization- forecast  run  was  made.   The  wind  data  were  taken 
from  the  data  at  12  hours  into  the  free  run.   Heights  in  the  tropics  were  re- 
placed by  a  constant  value  (admittedly  near  the  mean  of  the  "real"  data)  and  a 
"diabatic"  type  of  initialization  was  done  using  24  complete  forward  backward 
cycles  (as  in  the  case  of  Figures  7  and  7a).   For  this  run  the  results  are 
displayed  in  Figure  11,  Figures  12a, b,c,  and  Figures  13a, b,c.   Figure  11  is 
analogous  to  Figure  8,  Figures  9a,b,c  are  analogous  to  Figures  12a, b,c  and 
Figures  10a, b,c  are  analogous  to  Figures  13a,b,c.   These  results  clearly  show 
that  the  recovery  of  the  mass  field  by  the  initialization  and  the  subsequent 
forecast  are  quite  remarkable.   Note  that  some  information  regarding  the 
tropical  mass  field  was  needed,  if  only  just  its  mean  value.   The  procedure 
can  work  with  widely  divergent  values  of  the  original  (to  the  process)  mass, 
but  the  results  will  not  be  very  useful. 

The  remainder  of  the  figures  (14  through  25a, b,c  show  results  from  a  group 
of  runs  in  which  Rossby-gravity  type  waves  were  involved.   In  these  cases  all 
or  some  of  the  initial  fields  were  specified  analytically  and  no  mass  source 
or  friction  terms  were  used.   Specifically,  for  the  first  runs  for  which 
results  are.  given  in  Figures  14,  15a, b,c,  and  16a, b,c  the  initial  u,  v  and  $ 
fields  were  specified  from, 

v  =  A  sin  uX  e  s  (10) 

-r2/2 
u  =  Ky  A  cos  uX  e  '  (11) 


2 

4>  =  -  [—  +  v  Ky]  A  cos  uX  e  ^  '"  +  *  (12) 


37  ,  v  „_,  .  ___   .  -zTll 


9  /     1  //i  9      9 

where  u  -  ~,   L  =  18  AX,  A  =  10  m/sec,  ?  =  /By/ 4    ,  $  =  9.81  x  5400  m"/sec", 

2  /  2 

3  =  ~,   K  =  (B  -  *  $  C^)/^     i   -  v) ,  Cx  =  C/y,  v  =  *lA  /3(-  |  +  /\     +1)  and 

2  7T   1/4 

k  =  y~   ^3  $    .   Here  y=0  at  the  equator  and  is  +  the  distance  in  meters  from 
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Figure  11.   Forecast  run  with  forced  Kelvin  wave,  temporal  evolution  of  *,  U, 
V  at  selected  points  near  each  other  in  C  schema  at  2  S  (V  on 
equator)  ,  "diabatic"  dynamic  initialization. 


the  equator  respectively.   These  relations  either  come  from  or  can  be  derived 
from  those  given  in  Haltiner  and  Williams,  page  100,  following  the  original 
work  of  Matsuno  (1966).   Figure  14  again  displays  the  temporal  evolution  of  u, 
v,  $  at  the  point  at  2   degrees  South  and  Figures  15a, b,c  and  16a, b,c  are  the 
field  maps  at  hours  0  and  36  respectively. 

Next,  a  run  was  made  in  which  only  the  wind  fields  were  initially  specified 
from  Eqs  (10)  and  (11),  but  in  which  the  $  field  would  be  generated  by  the 
adiabatic  dynamical  initialization  procedure.   Then  a  36  hour  forecast  could 
be  made  using  these  fields.   Figures  17,  18a, b,c  and  19a, b,c  show  the 
analogous  data  for  this  run  to  that  shown  in  Figures  14  to  16.   Again  the  ini- 
tialization recaptures  the  mass  field  quite  well.   In  fact  the  total  amplitude 
of  the  mass  wave  computed  from  Eq  (12)  (Figure  15c)  is  slightly  greater  than 
that  from  the  initialization  (Figure  17c). 

Finally,  a  forecast  type  of  experiment  was  done  for  the  Rossby-gravity  type 
wave  as  was  done  before  for  the  Kelvin  wave  situation  and  the  earlier  cases 
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Figure  14.   Free  run,  Rossby-gravity  wave,  temporal  evolution  of  $ ,  U,  V. 
Fields  initialized  analytically. 


So  the  model,  using  all  the  Eqs  (10)-(12)  to  start  the  run,  was  integrated  for 
84  hours.   The  data  at  hour  48  were  stored  so  a  reinitialization  of  the  i 
field  could  be  done  at  that  point  and  a  36  hour  forecast  could  be  made  from 
that  data.   Figures  20,  21a, b,c  and  22a, b,c  show  results  from  the  free  run 
with  the  field  maps  shown  for  hour  48  and  hour  84.   Finally  Figures  23, 
24a, b,c  and  25a, b,c  show  the  results  for  the  forecast  run  with  the  field  maps 
shown  at  the  termination  of  the  initialization  and  at  36  hours. 

It  is  of  special  interest  to  compare  the  t>    field  in  Figure  21c  and  the 
<t>  fields  shown  in  Figure  24c.  The  field  in  Figure  21c  is  a  "snapshot"  of  the  <j> 
field  at  hour  48  in  the  regular  free  run.   To  initialize  the  "forecast"  run 
the  u  and  v  fields  at  hour  48  (Figures  21a,  21b)  were  retained,  a  constant 
value  for  $  from  30  N  to  30  S  was  prescribed,  and  the  initialization  procedure 
was  started.   By  its  termination  the  $  field  in  Figure  24c  was  generated.   The 
highest  and  lowest  values  are  recaptured  to  within  about  15  meters.   The  36 
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Figure  17.   Free  run,  Rossby-gravity  wave,  temporal  evolution  of  $,  U,  V. 
Initial  winds  obtained  analytically.   Initial  *'s  from  dynamic 
initialization. 


hour  forecast  $  field  (Figure  25c)  is  then  to  be  compared  to  the  84  hour  <j> 
value  in  the  control  or  free  run  (Figure  22c).   Now  the  high  and  low  centers 
are  within  10  meters.   Thus  this  method  of  initialization  shows  promise  in 
capturing  a  difficult  class  of  wave  motions  which  may  be  of  great  importance 
in  the  atmosphere  and  which  cannot  be  recaptured  easily  by  any  balancing  type 
procedure  -  at  least  right  near  the  equator.   These  waves  may  have  an  effect 
on  NWP  models  namely;  the  need  to  observe  and  initialize  them  properly  may  be 
absolutely  imperative  in  order  to  increase  large  scale  predictability  in  the 
atmosphere  beyond  4  or  5  days. 
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Figure  20.   Same  as  Figure  14  but  for  long  (84  hour)  free  run, 
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Figure  23.   Forecast  run,  Rossby-gravity  wave,  temporal  evolution  of  <j>,  U,  V. 
Initial  winds  from  hour  48  of  free  run.   Initial  <J>'s  obtained  from 
dynamic  initialization. 
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III.   Discussion 

A  dynamical  method  of  initialization  has  been  used  to  study  the  problem  of 
initializing  inertia-gravity  motions.   This  problem  is  not  of  such  high 
priority  at  the  present  time  in  NWP,  but  will  become  more  so  in  the  future 
because  interest  is  increasing  in  the  prediction  of  a)  meso-scale  phenomena  in 
midlatitudes  and  b)  tropical  motions  of  all  types.   In  the  latter  case  Kelvin 
and  Rossby-gravity  type  waves  at  or  near  the  equator  cause  problems  in  all 
methods  of  conventional  balance  type  initialization  and  objective  analyses 
schemes.   Here  it  is  demonstrated  chat  for  Kelvin  and  Rossby-gravity  waves  the 
dynamic  initialization  does  a  very  adequate  job  of  recapturing  the  mass  field 
near  the  equator.   For  operational  needs  however,  several  further  considera- 
tions would  have  to  be  investigated.   In  these  experiments  perfect  wind  obser- 
vations were  always  assumed.   To  relax  this  assumption  leads  into  the  area  of 
problems  related  to  objective  analyses-   In  this  connection  Anthes  (1976)  has 
proposed  in  his  work  on  hurricane  simulation  a  method  -  called  "nudging"  - 
useful  for  both  objective  analysis  and  initialization.   In  essence  this  method 
adds  another  coefficient  to  the  last  term  in  equation  (4)  so  it  would  be  of 
the  form,  C(a,x,y,t)  K   ( a^    -  a  )  where  the  C's  would  vary  from  grid  point 
to  grid  point  and  would  reflect  a  "confidence"  factor  in  the  observations. 
This  idea  could  be  incorporated  into  the  procedure  used  in  this  paper.   Also 
it  must  be  acknowledged  that  this  procedure  might  best  be  applied  in  conjunc- 
tion -  not  opposition  to  -  the  best  metnod  of  obtaining  balanced  winds  and 
mass  fields.   Certainly  this  would  be  of  great  aid  in  those  areas  between 
observed  data  points. 

Before  final  comments  a  brief  point  regarding  the  use  of  the  Euler  back- 
ward scheme.   If  a  standard  Von  Neumann  stability  analysis  is  carried  out  for 

9    /     2~~   4 
the  Euler  backward  scheme,  the  results  shows   X  '  =  /L-o  +  a     <  I  where 
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_  At   .   2ttAX 

°  =  c  a  sin  ~T~  ' 

In  the  situation  of  the  later  runs,  the  approximate  speed  of  the  pure  gravity 

wave,  Cg,  is  about  230  m/sec  with  At  =  360  sec  and  AX  =  444  kilometers.   There-1 

: 

fore  Cg(At/AX)  is  only  about  .2  near  the  equator.   Even  where  sin(2irAX/L)  =  1 

2    4  ~ 
(for  L  -   4  AX),  /l  -  o  +   a     =  .98  so  there  is  not  much  damping  in  these  runs 

for  these  waves.   For  smaller  values  of  C  or  longer  waves  up  to  L  =  18  AX  the 
damping  is  even  less  significant.   This  is  mentioned  to  emphasize  that  runs 
could  be  easily  designed  to  take  greater  advantage  of  the  damping  than  was 
done  for  the  cases  run  here. 

Finally,  it  may  be  of  interest  to  investigate  use  of  this  initialization 
procedure  using  other  types  of  models  such  as  spectral  or  finite  element 
models.   Though  there  is  no  physical  reason  to  believe  it  would  not  work,  these 
types  of  models  often  make  use  of  a  semi- implicit  scheme  which  slows  down  the 
inertia  gravity  waves  considerably.   How  this  might  affect  the  process  may  be 
of  interest.   Use  of  split-explicit  schemes  however  should  be  no  formal  prob- 
lem.  Here  the  external  gravity  mode  is  still  explicitly  predicted.   In  the 
case  of  the  spectral  formulation  restoration  coefficients  might  be  formed 
directly  on  the  basis  of  scale  considerations  relative  to  the  Rossby  radius  of 
deformation  in  both  directions.   All  these  areas  of  further  research  seem 
promising  . 
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